Intro

Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod

# Load libraries, install if needed
library(tidyverse)
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(qwraps2) 
library(wesanderson)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
library(brms)
library(sdmTMBextra)

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/density_model_cache/html")

For maps

# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
sf::sf_use_s2(FALSE)
# https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data

swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

# Define plotting theme for main plot
theme_plot <- function(base_size = 11, base_family = "") {
  theme_light(base_size = base_size, base_family = "") +
    theme(
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 9, colour = 'gray10', margin = margin(b = 1, t = 1)),
      strip.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 11, base_family = "") {
  theme_light(base_size = base_size, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 7),
        strip.text = element_text(size = 8, colour = 'gray10', margin = margin(b = 1, t = 1)),
        strip.background = element_rect(fill = "gray95"),
        legend.direction = "horizontal",
        legend.margin = margin(1, 1, 1, 1),
        legend.box.margin = margin(0, 0, 0, 0),
        legend.key.height = unit(0.4, "line"),
        legend.key.width = unit(2, "line"),
        legend.spacing.x = unit(0.1, 'cm'),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())
}

# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2

plot_map_raster <- 
ggplot(swe_coast_proj) + 
  xlim(xmin2, xmax2) +
  ylim(ymin2, ymax2) +
  labs(x = "Longitude", y = "Latitude") +
  geom_sf(size = 0.3) + 
  theme_plot() +
  guides(colour = guide_colorbar(title.position = "top", title.hjust = 0.5),
         fill = guide_colorbar(title.position = "top", title.hjust = 0.5))

Read data

d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/catch_q_1_4.csv")

# Calculate standardized variables
d <- d %>% 
  filter(year < 2020 & year > 1992 & quarter == 4) %>% 
  mutate(oxy_sc = oxy,
         temp_sc = temp,
         depth_sc = depth,
         ) %>%
  mutate_at(c("oxy_sc", "temp_sc", "depth_sc"),
            ~(scale(.) %>% as.vector)) %>% 
  mutate(year = as.integer(year)) %>% 
  drop_na(oxy, depth, temp)

plot_map_raster +
  geom_point(data = d, aes(X*1000, Y*1000))

Read the prediction grids

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")

pred_grid <- bind_rows(pred_grid1, pred_grid)

# Standardize data with respect to data grid:
pred_grid <- pred_grid %>%
  mutate(year = as.integer(year)) %>% 
  filter(year %in% c(unique(d$year))) %>% 
  mutate(depth_sc = (depth - mean(d$depth))/sd(d$depth),
         temp_sc = (temp - mean(d$temp))/sd(d$temp),
         oxy_sc = (oxy - mean(d$oxy))/sd(d$oxy)) %>% # Need to scale these to the mean and sd in the data!
  drop_na(oxy, depth, temp)

Make spde mesh

spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 200, 
                  type = "kmeans", seed = 42)

# Plot and save spde
png(file = "figures/supp/density/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()

Fit the density model

# Depth spline + oxy spline
# Takes about 30 minutes
m <- sdmTMB(density ~ 0 + as.factor(year) + s(depth_sc) + s(oxy_sc) + s(temp_sc),
            data = d,
            mesh = spde,
            family = tweedie(link = "log"),
            spatiotemporal = "AR1",
            spatial = "on",
            time = "year",
            reml = FALSE,
            control = sdmTMBcontrol(newton_steps = 1))
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.5.4
#> Current Matrix version is 1.5.3
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
#> glmmTMB was built with TMB version 1.9.3
#> Current TMB version is 1.9.4
#> Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)

tidy(m, conf.int = TRUE)

sanity(m)

Check QQ plot and residuals

samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
#> Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems

mcmc_res <- residuals(m, type = "mle-mcmc", mcmc_samples = samps)

png(file = "figures/supp/density/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()

d$residuals <- mcmc_res
# Residuals on map
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals), size = 0.5) +
  scale_colour_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  labs(color = "residuals") +
  theme_facet_map() +
  geom_sf(size = 0.1)
#> Warning: Removed 8 rows containing missing values (`geom_point()`).

  
ggsave("figures/supp/density/spatial_resid.png", width = 6.5, height = 8.5, dpi = 600)
#> Warning: Removed 8 rows containing missing values (`geom_point()`).

Check the AR1 parameter (rho is ar_phi on the -1 to 1 scale):

tidy(m, effects = "ran_pars", conf.int = TRUE)
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.5.4
#> Current Matrix version is 1.5.3
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> Standard errors intentionally omitted because they have been calculated in log
#> space.

Predict on grid

predict_mcod <- predict(m, newdata = pred_grid)

# Save prediction as csv (smaller than the model object...)
write.csv(predict_mcod, "output/predict_mcod_density.csv")

Plot predictions on map

# Plot predicted density and random effects
plot_map_raster +
  geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~ year, ncol = 5) +
  labs(fill = expression(kg/km^2)) +
  theme_facet_map() +
  geom_sf(size = 0.1)
#> filter: removed 42,282 rows (16%), 216,297 rows remaining

 
ggsave("figures/supp/density/est_map.png", width = 6.5, height = 8.5, dpi = 600)

# Plot spatiotemporal random effect
plot_map_raster +
  geom_raster(data = predict_mcod, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/density/epsilon_st_map.png", width = 6.5, height = 8.5, dpi = 600)

# Plot spatial random effect
plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
  scale_fill_gradient2() +
  geom_sf(size = 0.1)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining

  
ggsave("figures/supp/density/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)

Annual indices using get_index_sims

div_index_list <- list()

for(i in unique(pred_grid$sub_div)){
  
  pred_grid_sub <- pred_grid %>% filter(sub_div == i & depth > 10 & depth < 110)
  
  pred_sim_sub <- predict(m, newdata = pred_grid_sub, nsim = 500)
  
  index_sim <- get_index_sims(pred_sim_sub, 
                              area = rep(4*4, nrow(pred_sim_sub))) %>% 
    mutate(sub_div = i)
           
  div_index_list[[i]] <- index_sim
  
}
#> filter: removed 221,103 rows (86%), 37,476 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#> 
#> filter: removed 188,109 rows (73%), 70,470 rows remaining
#> 
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#> 
#> filter: removed 204,687 rows (79%), 53,892 rows remaining
#> 
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#> 
#> filter: removed 240,516 rows (93%), 18,063 rows remaining
#> 
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#> 
#> filter: removed 230,877 rows (89%), 27,702 rows remaining
#> 
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA

div_index <- bind_rows(div_index_list)

# In addition to the subdivision indicies, we'll also group a few and also do the total
# Now get the total index by specifying the area of a grid cell
# Total prediction 
preds_cod_sim <- predict(m, newdata = pred_grid, nsim = 500)

# Now calculate the biomass index
index_sim <- get_index_sims(preds_cod_sim, area = rep(4*4, nrow(preds_cod_sim)))
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.

# Join indices
index_sim <- index_sim %>% mutate(sub_div = "Total")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA

div_index_sim <- bind_rows(div_index %>% mutate(sub_div = as.character(sub_div)),
                           index_sim %>% mutate(sub_div = "Total"))
#> mutate: converted 'sub_div' from double to character (0 new NA)
#> mutate: no changes

Compare the above biomass-index with one where I filter deep depths

# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 10 & depth < 120) %>% nrow()
#> filter: removed 250,568 rows (97%), 8,011 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 10 & depth < 120)
#> filter: removed 42,282 rows (16%), 216,297 rows remaining

pred_avg_sim_sub <- predict(m, newdata = pred_grid_sub, nsim = 500)

index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(4*4, nrow(pred_avg_sim_sub)))
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.

# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_sim, area = "full"),
                                mutate(index_avg_sim_sub, area = "subset")) %>% 
  mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'est_t' (double) with 54 unique values and 0% NA
#>         new variable 'lwr_t' (double) with 54 unique values and 0% NA
#>         new variable 'upr_t' (double) with 54 unique values and 0% NA

ggplot(index_avg_sim_comp, aes(year, est_t, ymin = lwr_t, ymax = upr_t, color = area)) +
  geom_line() +
  geom_ribbon(alpha = 0.4)

Combined plot with biomass index by sd, total and predictions on map for two years

lab_size <- 1.7

div_index_sim <- div_index_sim %>% mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) 

ind_plot_sd <- 
  ggplot() +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  geom_line(data = filter(div_index_sim, !sub_div == "Total"),
            aes(year, est_t, color = sub_div)) + 
  geom_ribbon(data = filter(div_index_sim, !sub_div == "Total"),
              aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes") +
  theme_plot() + 
  guides(colour = "none",
         fill = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme(legend.position = c(0.6, 0.9),
        legend.direction = "horizontal",
        legend.key.height = unit(0.03, "cm"),
        legend.key.width = unit(0.4, "cm"),
        legend.text = element_text(size = 6),
        legend.title = element_text(size = 7),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        panel.background = element_blank())

ind_plot_tot <- 
  ggplot() +
  geom_line(data = filter(div_index_sim, sub_div == "Total"),
            aes(year, est_t)) +
  geom_ribbon(data = filter(div_index_sim, sub_div == "Total"),
              aes(year, ymin = lwr_t, ymax = upr_t, fill = "Total"), alpha = 0.4) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_fill_manual(values = "grey40") + 
  labs(x = "Year", y = "Tonnes", fill = "") +
  theme_plot() + 
  theme(legend.position = c(0.88, 0.92),
        legend.key.height = unit(0.03, "cm"),
        legend.key.width = unit(0.4, "cm"),
        legend.text = element_text(size = 6),
        legend.title = element_text(size = 7),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

cpue_map_94_18 <- plot_map_raster +
  geom_raster(data = filter(predict_mcod, year %in% c("1994", "2018") & depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~ year, ncol = 2) +
  labs(fill = expression(kg/km^2)) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.5, "cm"),
        legend.position = c(0.9, .085),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA),
        legend.text = element_text(size = 5),
        legend.title = element_text(size = 7)) +
  guides(fill = guide_colorbar(title.position = "top", title.hjust = 0.9)) +
  geom_sf(size = 0.1) +
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)

(ind_plot_sd | ind_plot_tot) / cpue_map_94_18 + plot_annotation(tag_levels = 'A')


ggsave("figures/density.pdf", width = 17, height = 15, units = "cm")

# All years cpue
plot_map_raster +
  geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~ year, ncol = 5) +
  labs(fill = expression(kg/km^2)) +
  theme_facet_map() +
  geom_sf(size = 0.1)


ggsave("figures/supp/density/all_years_density_pred.png", width = 6.5, height = 8.5, dpi = 600)

Biomass weighted quantiles of depth, oxygen and temperature

# Plot bathymetry
plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = depth)) +
  scale_fill_viridis(option = "A", direction = -1) +
  labs(fill = "Depth") +
  theme(plot.margin = unit(c(0, 0, 0.5, 0), "cm"),
        legend.key.height = unit(0.1, "cm"),
        legend.key.width = unit(0.7, "cm"),
        legend.position = c(0.8, .07),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA)) +
  geom_sf(size = 0.1) +
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining

lab_size <- 1.6
## Depth
# Plot bathymetry
bathy_plot <- plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = depth)) +
  scale_fill_viridis(option = "A", direction = -1) +
  labs(fill = "Depth", x = "", y = "") +
  theme_plot(base_size = 11) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.key.height = unit(0.1, "cm"),
        legend.key.width = unit(0.35, "cm"),
        legend.position = c(0.8, .085),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA),
        legend.text = element_text(size = 5),
        legend.title = element_text(size = 6)) +
  guides(fill = guide_colorbar(title.position = "top", title.hjust = 1)) +
  geom_sf(size = 0.1) +
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining
  
# Now calculate quantiles of the biomass-weighted values
pal <- c(rev(brewer.pal(n = 8, name = "Dark2")), "gray20")

wm_depth <- predict_mcod %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)),
            "1st" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.25)),
            "3rd" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.75))) %>% 
  pivot_longer(cols = c("Median", "1st", "3rd"),
               names_to = "series", values_to = "depth") %>% 
  group_by(series) %>%  # standardize within for easy plotting
  mutate(depth_ct = depth - mean(depth)) %>% 
  ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st, 3rd) into (series, depth) [was 27x4, now 81x3]
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'depth_ct' (double) with 79 unique values and 0% NA
#> ungroup: no grouping variables
  
plot_w_dep <- ggplot(wm_depth, aes(year, depth, color = series, group = series, fill = series)) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 0.8) +
  geom_point(size = 1.3, alpha = 0.8) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_color_manual(values = pal, name = "Quantiles") +
  scale_fill_manual(values = pal) +
  guides(fill = "none", color = "none") +
  labs(y = "Depth [m]", x = "Year", color = "") +
  scale_y_reverse() + 
  theme_plot(base_size = 11) +
  theme(plot.margin = unit(c(0.9, 0, 0, 0), "cm"))
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

## Oxygen
# Plot oxygen in 2006
oxy_plot <- plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = oxy)) +
  scale_fill_viridis() +
  labs(fill = expression(paste("O" [2], " [ml/L]", sep = "")), x = "") +
  theme_plot(base_size = 11) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.key.height = unit(0.1, "cm"),
        legend.key.width = unit(0.35, "cm"),
        legend.position = c(0.8, .085),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA),
        legend.text = element_text(size = 5),
        legend.title = element_text(size = 6)) +
  guides(fill = guide_colorbar(title.position = "top", title.hjust = 1)) +
  geom_sf(size = 0.1) +
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining

# Plot biomass-weighted oxygen
# Weighted total
wm_oxy <- predict_mcod %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
            "1st" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
            "3rd" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>% 
  pivot_longer(cols = c("Median", "1st", "3rd"),
               names_to = "series", values_to = "oxy") %>% 
  ungroup() %>% 
  group_by(series) %>% # standardize within for easy plotting
  mutate(oxy_ct = oxy - mean(oxy)) %>%
  ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st, 3rd) into (series, oxy) [was 27x4, now 81x3]
#> ungroup: no grouping variables
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'oxy_ct' (double) with 81 unique values and 0% NA
#> ungroup: no grouping variables

# Find which depths to calculate average oxygen on
wm_depth %>% filter(series == "1st") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
wm_depth %>% filter(series == "3rd") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped

annual_oxy <- pred_grid %>%
  drop_na(oxy) %>%
  filter(depth > 29 & depth < 61) %>% ###
  #filter(depth > 40 & depth < 60) %>% ###
  group_by(year) %>%
  summarize(median_oxy = median(oxy)) %>%
  ungroup() %>% 
  rename("oxy" = "median_oxy") %>% 
  mutate(series = "Median (env.)",
         oxy_ct = oxy - mean(oxy))
#> drop_na: no rows removed
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> rename: renamed one variable (oxy)
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'oxy_ct' (double) with 27 unique values and 0% NA
  
oxygen_series <- bind_rows(wm_oxy, annual_oxy)

plot_w_oxy <-
ggplot(oxygen_series, aes(year, oxy, color = series, group = series, fill = series, linetype = series)) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 0.8) +
  geom_point(size = 1.3, alpha = 0.8) + 
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_linetype_manual(values = c(1, 1, 1, 2)) +
  scale_color_manual(values = pal, name = "Quantiles") +
  scale_fill_manual(values = pal) +
  guides(fill = "none", linetype = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0)) +
  labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
       color = "") +
  theme_plot(base_size = 11) +
  theme(legend.position = c(0.25, 0.12),
        legend.key.size = unit(0.25, "cm"),
        legend.background = element_rect(fill = NA),
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

plot_w_oxy


# zoom in..
plot_w_oxy + 
  scale_y_continuous(breaks = seq(0, 10, by = 0.1))

  #coord_cartesian(ylim = c(4.5, 5))
#

oxygen_series %>% as.data.frame()

## Temperature
# Plot temperature in 2006
temp_plot <- plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = temp)) +
  scale_fill_viridis(option = "B", direction = -1) +
  labs(fill = "Temperature [°C]", y = "") +
  theme_plot(base_size = 11) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.key.height = unit(0.1, "cm"),
        legend.key.width = unit(0.35, "cm"),
        legend.position = c(0.8, .085),
        legend.direction = "horizontal",
        legend.background = element_rect(fill = NA),
        legend.text = element_text(size = 5),
        legend.title = element_text(size = 6)) +
  guides(fill = guide_colorbar(title.position = "top", title.hjust = 1)) +
  geom_sf(size = 0.1) +
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining

# Plot biomass-weighted temperature
# Weighted total
wm_temp <- predict_mcod %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.5)),
            "1st" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.25)),
            "3rd" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.75))) %>% 
  pivot_longer(cols = c("Median", "1st", "3rd"),
               names_to = "series", values_to = "temp") %>% 
  ungroup() %>% 
  group_by(series) %>% # standardize within for easy plotting
  mutate(temp_ct = temp - mean(temp)) %>%
  ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st, 3rd) into (series, temp) [was 27x4, now 81x3]
#> ungroup: no grouping variables
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'temp_ct' (double) with 81 unique values and 0% NA
#> ungroup: no grouping variables

# Find which depths to calculate average oxygen on
wm_depth %>% filter(series == "1st") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
wm_depth %>% filter(series == "3rd") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped

annual_temp <- pred_grid %>%
  drop_na(temp) %>%
  filter(depth > 29 & depth < 61) %>% ###
  #filter(depth > 40 & depth < 60) %>% ###
  group_by(year) %>%
  summarize(median_temp = median(temp)) %>%
  ungroup() %>% 
  rename("temp" = "median_temp") %>% 
  mutate(series = "Median (env.)",
         temp_ct = temp - mean(temp))
#> drop_na: no rows removed
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> rename: renamed one variable (temp)
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'temp_ct' (double) with 27 unique values and 0% NA
  
temp_series <- bind_rows(wm_temp, annual_temp)

plot_w_temp <-
ggplot(temp_series, aes(year, temp, color = series, group = series, fill = series, linetype = series)) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 0.8) +
  geom_point(size = 1.3, alpha = 0.8) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_linetype_manual(values = c(1,1,1,2)) +
  scale_color_manual(values = pal, name = "Quantiles") +
  scale_fill_manual(values = pal) +
  guides(fill = "none", linetype = "none", color = "none") +
  labs(y = "Temperature [°C]", x = "Year", color = "") +
  theme_plot(base_size = 11) +
  theme(plot.margin = unit(c(0, 0, 0.9, 0), "cm"))

plot_w_temp


## Make combined plot

(bathy_plot | oxy_plot | temp_plot) / (plot_w_dep | plot_w_oxy | plot_w_temp) +
  plot_annotation(tag_levels = 'A') &
  theme(axis.text = element_text(size = 7),
        axis.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        legend.title = element_text(size = 6))


l <- (bathy_plot / oxy_plot / temp_plot) & theme(legend.text = element_text(size = 6),
                                                 legend.title = element_text(size = 7))

r <- (plot_w_dep / plot_w_oxy / plot_w_temp) & theme(legend.text = element_text(size = 8),
                                                     legend.title = element_text(size = 9))

(l | r) + plot_annotation(tag_levels = 'A') 

  

#ggsave("figures/weighted_quantiles.pdf", width = 17, height = 11, units = "cm")
ggsave("figures/weighted_quantiles.pdf", width = 17, height = 20, units = "cm")

Weighted oxygen and depth by year and sub div

# Oxygen
wm_oxy_sd <- predict_mcod %>%
  filter(depth > 29 & depth < 61) %>%
  group_by(year, sub_div) %>%
  summarise("weighted median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
            "un-weighted median" = median(oxy)) %>% 
  pivot_longer(c("weighted median", "un-weighted median"))
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 4 columns, one group variable remaining (year)
#> pivot_longer: reorganized (weighted median, un-weighted median) into (name, value) [was 135x4, now 270x4]

ggplot(wm_oxy_sd, aes(year, value, color = factor(name), fill = factor(name))) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 0.8) +
  geom_point(size = 1.5) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_color_brewer(palette = "Accent", name = "") +
  scale_fill_brewer(palette = "Accent") +
  guides(fill = "none") +
  labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year") +
  theme_plot() +
  facet_wrap(~sub_div) +
  theme(aspect.ratio = 1,
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0.2), "cm"))


ggsave("figures/supp/density/weighted_oxygen_subdiv.png", width = 6.5, height = 6.5, dpi = 600)

# Calculate weighted mean and quantiles for sd 25 only 
predict_mcod %>%
  filter(sub_div == 25 & year >= 2015) %>%
  summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>%
  as.data.frame()
#> filter: removed 245,134 rows (95%), 13,445 rows remaining
#> summarise: now one row and 3 columns, ungrouped

# Now compare delta change in oxygen with delta change in condition

# Depth
wm_dep_sd <- predict_mcod %>%
  filter(depth > 29 & depth < 61) %>%
  group_by(year, sub_div) %>%
  summarise("weighted median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)),
            "un-weighted median" = median(depth)) %>% 
  pivot_longer(c("weighted median", "un-weighted median"))
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 4 columns, one group variable remaining (year)
#> pivot_longer: reorganized (weighted median, un-weighted median) into (name, value) [was 135x4, now 270x4]

ggplot(wm_dep_sd, aes(year, value, color = factor(name), fill = factor(name))) +
  stat_smooth(data = filter(wm_dep_sd, name == "weighted median"), method = "gam", formula = y ~ s(x, k = 4), size = 0.8) +
  geom_point(size = 1.5) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_color_brewer(palette = "Accent", name = "") +
  scale_fill_brewer(palette = "Accent") +
  guides(fill = "none") +
  labs(y = "Depth [m]", x = "Year", color = "") +
  theme_plot() +
  facet_wrap(~sub_div) +
  theme(aspect.ratio = 1,
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
#> filter (grouped): removed 135 rows (50%), 135 rows remaining


ggsave("figures/supp/density/weighted_depth_subdiv.png", width = 6.5, height = 6.5, dpi = 600)

# Now plot together
p1 <- ggplot(filter(wm_dep_sd, name == "weighted median"),
       aes(year, value, color = factor(sub_div))) +
  stat_smooth(data = filter(wm_dep_sd, name == "weighted median"), method = "gam", formula = y ~ s(x, k = 4), size = 0.8, se = FALSE) +
  geom_point(size = 1.5) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_color_brewer(palette = "Dark2", name = "") +
  guides(fill = "none") +
  labs(y = "Depth [m]", x = "Year", color = "") +
  theme_plot() +
  theme(aspect.ratio = 1,
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
#> filter (grouped): removed 135 rows (50%), 135 rows remaining

wm_oxy_sd2 <- filter(wm_oxy_sd, name == "weighted median")
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
m1 <- lm(value ~ year*as.factor(sub_div), data = wm_oxy_sd2)
m2 <- lm(value ~ year+as.factor(sub_div), data = wm_oxy_sd2)

p2 <- ggplot(wm_oxy_sd2, aes(year, value, color = factor(sub_div))) +
  stat_smooth(data = filter(wm_oxy_sd, name == "weighted median"), method = "gam", formula = y ~ s(x, k = 4), size = 0.8, se = FALSE) +
  geom_point(size = 1.5) +
  scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
  scale_color_brewer(palette = "Dark2", name = "") +
  guides(fill = "none") +
  labs(y = "Depth [m]", x = "Year", color = "") +
  theme_plot() +
  theme(aspect.ratio = 1,
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
#> filter (grouped): removed 135 rows (50%), 135 rows remaining

p2

Conditional effect of oxygen

# Prepare prediction data frame
d2 <- d %>% drop_na(oxy)
nd_oxy <- data.frame(oxy = seq(min(d2$oxy), max(d2$oxy), length.out = 100))

nd_oxy <- nd_oxy %>%
  mutate(year = 2003L,
         depth_sc = 0,
         oxy_sc = (oxy - mean(oxy))/sd(oxy),
         temp_sc = 0)

# Predict
p_cond_oxy <- predict(m, newdata = nd_oxy, se_fit = TRUE, re_form = NA)

cond_oxy <- ggplot(p_cond_oxy, aes(oxy, exp(est),
  ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(30, 220)) +
  theme(aspect.ratio = 1) +
  labs(x = expression(Oxygen[haul]))

Conditional effect of depth

# Prepare prediction data frame
nd_dep <- data.frame(depth = seq(min(d2$depth), max(d2$depth), length.out = 100))

nd_dep <- nd_dep %>%
  mutate(year = 2003L,
         depth_sc = (depth - mean(depth))/sd(depth),
         oxy_sc = 0,
         temp_sc = 0)

# Predict
p_cond_dep <- predict(m, newdata = nd_dep, se_fit = TRUE, re_form = NA)

cond_depth <- ggplot(p_cond_dep, aes(depth, exp(est),
  ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(30, 220)) +
  theme(aspect.ratio = 1) +
  labs(x = expression(Depth[haul]))

Conditional effect of temperature

# Prepare prediction data frame
nd_temp <- data.frame(temp = seq(min(d$temp), max(d$temp), length.out = 100))

nd_temp <- nd_temp %>%
  mutate(year = 2003L,
         depth_sc = 0,
         oxy_sc = 0,
         temp_sc = (temp - mean(temp))/sd(temp))

# Predict
p_cond_temp <- predict(m, newdata = nd_temp, se_fit = TRUE, re_form = NA)

cond_temp <- ggplot(p_cond_temp, aes(temp, exp(est),
  ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(30, 220)) +
  theme(aspect.ratio = 1) +
  labs(x = expression(Temp[haul]))

Plot together

p_cond_dep2 <- p_cond_dep %>%
  mutate(variable = "Depth") %>% 
  dplyr::select(est, est_se, depth_sc, variable) %>% 
  rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_oxy2 <- p_cond_oxy %>%
  mutate(variable = "Oxygen") %>% 
  dplyr::select(est, est_se, oxy_sc, variable) %>% 
  rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_cond_tem2 <- p_cond_temp %>%
  mutate(variable = "Temperature") %>% 
  dplyr::select(est, est_se, temp_sc, variable) %>% 
  rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

conds <- bind_rows(p_cond_dep2, p_cond_oxy2, p_cond_tem2)

ggplot(conds, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
                  ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  facet_wrap(~variable, ncol = 2) +
  labs(y = "Biomass density",
       x = "Scaled value") +
  theme_plot()


ggsave("figures/supp/density/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()